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Abstract We discuss the nonlinear stability of phase-locked states for globally coupled nonlinear 
oscillators with finite inertia, namely the modified Kuramoto model, in the context of the robust 
£°°-norm. We show that some classes of phase-locked states are orbitally £°°-stable in the sense that 
its small perturbation asymptotically leads to only the phase shift of the phase-locked state from 
the original one without changing its fine structures as keeping the same suitable coupling strength 
among oscillators and the same natural frequencies. The phase shift is uniquely determined by the 
average of initial phases, the average of initial frequencies, and the strength of inertia. We numerically 
confirm the stability of the phase- locked state as well as its uniqueness and the phase shift, where 
various initial conditions are considered. Finally, we argue that some restricted conditions employed 
in the mathematical proof are not necessary, based on numerical simulation results. 

Keywords Kuramoto oscillators, inertia, phase-locked state, orbital stability, phase shift 



1 Introduction 

Since Kuramoto introduced a mathematical model of coupled nonlinear oscillators [T] by refining the 
earlier Winfree's model j2] to be more mathematically tractable, it has become a minimal model for 
collective synchronization phenomena, which are ubiquitous in real systems ranging from physics to 
biology. The original Kuramoto model is very simple but exhibits lots of rich behaviors including 
a synchronization transition, where all the oscillators' phases are tuned by the coupling strength 
against the diversity of natural frequencies, and eventually reach a phase-locked state (frequency 
entrainment) including in-phase synchronization with exactly the same value (see [31[3] for detailed 
discussion) . 
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Synchronization is a quite interesting nonlinear phenomenon, which can yield a dynamic phase 
transition against the coupling strength among weakly coupled oscillators for a given natural fre- 
quency distribution. Above the coupling strength threshold, oscillators become partially phase- 
locked, which is called partial synchronization. The existence of partially phase-locked solutions 
and their stability were studied by Kuramoto [T] and Crawford [S], respectively. As an extension of 
the aforementioned works, Aeyels and Rogge addressed the stability issue for networks of a finite 
number of oscillators [6] , and MiroUo and Strogatz analyzed the linear stability of phase- locked states 
for globally coupled Kuramoto oscillators [3- Later on, lots of phase models have been proposed to 
describe the dynamic behavior of large populations of nonlinearly coupled oscillators. Furthermore, 
it has also been widely discussed that the nature of synchronization transitions for the continuum 
version of the original Kuramoto model is subject to the shape of the natural frequency distribution 
function. While the synchronization transition has been widely discussed from the physical point of 
view, the strict shape of phase-locked states has not drawn much less attention because it concerns 
far from the transition. Of course, the transition is one of the interesting topics to be studied in 
many senses, such as critical behaviors and its universality issue. However, it is mostly restricted to 
numerical studies, which implies that there are always some technical difficulties to reach a correct 
conclusion without rigorous mathematical guidelines and/or proofs. In that sense, it is meaningful 
if we can report what and how can be mathematically proven in the perspective of relevant aspects 
directly and indirectly to major interests for a given model before proceeding numerical studies near 
and at the transition. 

In this paper, we address the modified Kuramoto model with finite inertia and take an interest 
with the details of phase-locked states in order to handle the shape of phase-locked states as well as 
their uniqueness issue for a fixed distribution of natural frequencies. Consider an ensemble of finite 
number of Kuramoto oscillators with finite inertia. Let 9i = 8i {t) be the phase of the z-th oscillator 
and m denotes the strength of uniform inertia acting on Kuramoto oscillators. In this situation, the 
dynamics of 9i is governed by the initial valule problem of second-order ODEs: 

K ^ 

+ = r?, + — ^sin(0j-0,), i>0, i = l,...,iV, (1.1) 
subject to initial phase-frequency configuration: 

{e,,e,m = {e,o,Lo,o^ko)- (1.2) 

Without loss of generality, we assume that the average of natural frequencies is zero: 

1 ^ 

1=1 

The system (|l.ip was first introduced by Ermentrout [5] as a phenomenological model to explain the 
slow synchronization of certain biological systems, e.g., fircfiies of the Pteroptyx malaccae, and this 
model has been used to describe various dynamical systems [9l [T0l[m[T2ll 1 3 llMll 1 5 II 1 6| . The whole 
review and mathematical results for the Kuramoto model can be found in [^I fTTlfTSlIT^ . In the 
earlier study, the first author and his collaborators [ID] have verified that for some classes of initial 
configurations and mK > j, the speed to the phase- locked states is slower than that of the original 
Kuramoto model (without the inertial term), which illustrates the slow relaxation of Kuramoto 
oscillators. 

The main novelty of this paper is to provide a simple proof for the nonlinear stability of some 
classes of phase-locked states in globally coupled Kuramoto oscillators with finite inertia. Using 
the robust £°°-norm, we present the proof under the suitable conditions on initial phase-frequency 
configurations, coupling strengths, and finite inertia. It is noted that phase-locked states are orbitally 
£°°-stable, which results in the uniqueness of phase-locked states. 

This paper is organized after introduction as follows. In Sec. 2, we briefly recall some basic 
definitions and mathematical structures od the system (|l.ip . In Sec. 3, we show the proof of orbital 
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stability of phase-locked states and provide extensive numerical simulations, which confirm our 
analytical results. Finally, Section 4 is devoted to the summary of main results and the conclusion 
of the paper. 



2 Frameworks and main results 

In this section, we discuss two frameworks depending on the relative size of inertial strength and 
present main results on the orbital stability of phase-locked states. We also present several mathe- 
matical structures of the system (|l.ip - (|1.3p . and recall a second-order Gronwall's lemma from PO] 
which will be used in the proof of our main result. 

Before we begin our technical discussions, we first recall several definitions for the phase-locked 
state and its orbital stability as follows. 

Definition 2.1 (Phase-locked state and orbital stability) 

1. Lete{t) := (6*1 (t), • • • ,6'jv(t)) € be the solution to the system ([TTI) - ^T^. Then 9 is (strongly) 
phase-locked if and only if 9 is an equilibrium solution to the system of ODEs i.e., 



2. Let 9 and 9 be two phase-locked solutions to the system p.ip and (|1.3p . Then, 9 and 9 are 
congruent and denoted if and only if there exists a constant /3 G i? such that 



3. Let 9'^ be a phase-locked solution to the system (jl.ip and (jl.3p . Then 9^^ is orbitally stable if 
and only if for any initial phase- configuration 9q close to 9'^ in a norm \\ ■ \\, the solution 9{t) 
converges to the phase-shift of 9'' in a norm \ \ ■ \\ as t ^ oo, i.e.. 



Below we add some comments on the above definitions. 
Remark 2.1 

1. In general, a phase-locked state is a traveling profile with a constant phase velocity Qc- 

2. The time- dependent solution 9 = 9{t) is a (weakly) phase-locked state if and only if there exist 
positive constants C, , C* > independent of t satisfying 



9-9 = PIn 



where 1m (1, ■ ■ ■ ; 1) G 



lim \\9{t) - {9" -f ^Iat)!! = 0, for some constant (3. 



C,<\9,{t)-9j{t)\<C*, t>Q. 



2.1 Basic estimates 



We rewrite the system (jl.ip as a system of first-order ODEs for 



{9i, uji :— 9i): 



■I 1 



t>Q, i = l,--- ,N, 



1 



K 



N 



(2.1) 
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We next introduce macro variables (center-of-mass frame) and micro variables (fluctuations around 
macro variables): 



N , N 



1=1 1=1 

Then, macro- variables and micro- variables satisfy, respectively, 



g „ . _ uJc 

m 



(2.2) 



and 



K ^ 

§1 = LOt, niLJi = -uji + ■^J + ^ X] ^ ^'O- (2-3) 

By direct calculations, 

Gcii) = ^c(O) + mwe(0)(l - e-™), We(0 = e- ™ Wc(0). 

Thus, macro-variables converge toward some constant states that are determined by their initial 
configurations and the magnitude of inertia: 

lim (0,(<), wc(<)) = (^c(O) + mw,(0), 0). (2.4) 

For convenience, we recall the following second-order differential inequality: 

ay + fey + CT/ + d < 0, t > 0, 

y(0) = yo, 2/(0) -yi, ^ ' ^ 

where a > 0, 6, c and d are constants. 

Lemma 2.1 [20j Let y = y{t) be a nonnegative C'^ -function satisfying the differential inequality 
(|2.5p . Then we have following relations: 



{i) - 4ac > 0; 
(m) 6^ - 4ac < 0; 



, yi + i^iyo + 



6 , r 4a(i / 6 2d\ 

y(i)<e--*[yo + ^ + (^yo + 2/i + y)t 

where decay exponents vi and V2 are given as follows. 



Aad 



b + y/b'^ - Aac b - s/W^Aac 



■■= , 1^2 



2a 2a 
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2.2 Main results 

In order to present main results regarding the stability of phase-locked states, we first discuss two 
frameworks, which depend on the relative magnitude of inertia m to the strength of coupling K 
between oscillators. 

Let us recall a ^^'-norm for a finite-dimensional vector space. For 9 G T^, the £P-norm of 9 is 
defined as follows. 

N 1 



i=l 

max I I , p = oo. 



^ l<i<N ' 

Then, phase and natural frequency diameters can be denoted as 



5^i0.rr,pe [i,oo). 



D{9):^ max \9,-9j\ and D{Q) := max \fl,-nj\. 

It is noted that for 6 G with ^f^i Oi = 0, ||^||^^ and D{0) arc equivalent in the sense that 

\\9\\,^<D{9)<2\\9\\,^. 

The following two frameworks that concern the magnitude of inertia were first introduced in |20j 
for the existence of phase- locked states to the system (|l.ip . 

• Framework A: (Small inertia regime) 

1. The strength of coupling K and the magnitude of inertia m satisfy 



0< D{Q) < K, mK < 



4sinL»°° 

where £ ('^' ^) ^^'^ ^'^'^^ '^^ ^^^'^ following trigonometric equation: 

Dm) 

sm X = . 

K 

2. An initial configuration of {9q,ujo) satisfies 

<maxiD(9o),D(9o) + 2mD(9(t))\ ] <D°°, 
L \t=o) 

• Framework B: (Large inertia regime) 

1 . The strength of coupling K and the magnitude of inertia m satisfy 

2. An initial configuration of {9q,ujo) satisfies 

< maxlD(eo), D(9o) + 2mD(e(t)) ] < imDm). 
I t=aJ 
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< -. 



4 sin 



By definitions, Framework A and Framework B correspond some restricted cases of mK. Let V be 
the collection of all phase-loeked states evolved from various initial configurations satisfying either 
Framework A or Framework B, i.e., 

V {9" = {91, • • • , 61^) e : 9^ -.^ lim 9{t), where 9{t) is the solution to the system (HI]) - ([T3)) 

t—^oo 

with initial datum satisfying either Framework A or Framework B }. 

The set T' is a proper subset of all possible phase-locked states for the system (jl.ip and (jl.Sp (see 
Section 3.2), and phase-locked states in V have a diameter strictly less than ^. 

We now move onto the main results of this paper, the orbital stability of phase-locked states. 
Theorem 2.1 (Orbital stability) 

Suppose either Framework A or Framework B hold, and let 9'^ he a given phase-locked state in V . 
Then 9"^ is orbitally £°° -stable in the sense that for any perturbed initial configuration 9q satisfying 
the condition (2) in either Framework A or Framework B, the perturbed solution 9{t) satisfies 

lim \\9{t) - (r + (Z\0)°°ljv)||f- = 0, 

i— >oo 

where the phase-shift {A9)°^ is explicitly given by 

{A9r := 01 - OM - mCoM- 

Remark 2.2 

1. The conditions of Framework A and Framework B are independent of the system size N. Hence, 
our results can be lifted to the kinetic regime via the thermodynamic limit. 

2. Under both Framework A and Framework B, the initial phase- frequency configuration (9q,ujo) 
evolves toward the asymptotic phase-locked state (0°°,O); 

+ ^ E - = 0, sup D{9{t)) < f 

i—i — 

3. For the case of m = (the original Kuramoto model), the orbital stability of phase-locked states 
is provided in I21f using the -contraction theory: 

lim \\9it) ~ (r + {A9)°°In)\\ii = 0. 

t—^oo 

However, we cannot apply this estimate to the system (jl.ip due to inertia. Thus, we employ a new 
estimate to include the previous result and also cover the inertial effect, which is based on i°° -metric. 

4. Theorem 2.1 implies that a phase-locked state with the phase diameter D{9) .strictly less than ^ 
is unique up to the phase shift. Let 9 and 9 be the two phase-locked states emerged from initial data 
{9o,uio) and (0o,wo), respectively. Suppose that 

D{9°°), D{9'^) < |. 

Then, by the same argument as in Theorem 2.1, we have 

9^-9°^ = (^9,{0) - 9,{0) + m(we(0) - c:'c(0)))lw. 
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3 Nonlinear stability of phase-locked states 

In this section, wc claim the orbital nonlinear stability of the phase-locked states in the £°° norm as 
providing some analytical proof, which is also numerically confirmed. 



3.1 The proof of Theorem 2.1 

Suppose either Framework A or Framework B hold, and let 9*^ be a given phase-locked state in V 
with D{9^) < ^ and 9o be a perturbed phase-configuration of 9"^. Then, thanks to Theorem 2.1, the 
perturbed configuration 9 ~ 9{t) with the initial configuration 9q satisfies 



D{9'')+D{9{t)) KIT, t>0. 



We set 



N 



1 

Iv ^^^^^ ~ 



Am max cii, d™ min di, D(a(t)) ~ Am — Q^n 



Note that 



a,{t) = - = - 0,(0) - 771(1 - e--)a;,(0), 

^6';;-^c(0)-mWc(0), as t -> oo, (3.1) 
\aj(t) - a,{t)\ = \9] - 9j{t) - {91 - 9,{t))\ < D{9'') + D{9{t)) < tt, for t > 0. 

By simple calculations, we obtain 

+ ^ = ^ V cos + ^) sm (^) . (3.2) 



m- 



• Step A (Derivation of Gronwall's inequality for D{a)): 
It follows from that 

d'^&M , doiM 2K ^ (91 - 91 9j-9i\ . / d^ - d^/ \ 



df^ dt N 



^ / i^sin2i^g" \ ^ / d, -dM \ (3.3) 



A' sin 2Dg" 
2L>^'' 



Am, 



where we used the fact — tt < dj — dj\/ < to find 



i^(0^) + i^(0o) ^ 
1=1 

dj - dji/^ ^ 1^ sin D"^"" ^ dj - um ^ 



2 



Similarly, we find 



d^d„, dd„, A'sin2£>n'' 
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Table 1 All the conditions and parameter sets used in our simulation tests for random 100 realizations of 
9{0) G [-tt/IO, tt/IO] and ij{0) G [-Tr.Tr]: 





Distribution 


i'Vamcwork 


m 


K 


Figl,2(a)left 


uniform ( 


-1,1 


) 


A 


0.1 


2.0 


Figl,2(a)right 


uniform ( 


-1,1 


) 


B 


0.2 


2.0 


Figl,2(b)left 


Caucliy(7 — 1 


) 


A 


0.004 


40.0 


Figl,2(b)riglit 


Caucliy(7 — 1) 


B 


0.01 


40.0 


Fig4,5(a)left 


uniform ( 


-1,1 


) 




1.0 


3.0 


Fig4,5(a)right 


uniform( 


-1,1 


) 




0.2 


2.0 


Fig4,5(b)left 


Caucliy(7 — 1 


) 




0.5 


40.0 


Fig4,5(b)riglit 


Cauchy(7 — 1 


) 




0.1 


20.0 



Wc combine the estimates p.3p and p. 41) to find 



d^Dia) , dDia) , ^^^^^ < ^ ^sin2^^ ^3^^^ 



dt^ dt V ; - , • ^j^av 



• Step B (Decay estimates of D{a)): 

We apply Lemma 2.1 for (|3.5p to obtain 



D{a{t)) < 



D{ao)e-^'^* + 7'"' fi>(ao) + f^iD{ao)) , 1 - 4mA' > 0, 

\/l—4mK V / 



e 2„ 



D{ao) + (^^D{ao) + D{ao))t 



1 - AmK < 0, 



where 



_ 1 + Vl - irnK _ 1 - Vl - 4mK 



Hence, for any e e (o, 2m) ' have 

D(a(t)) « C'(l)e-^('')* for large timet, Afe) := min (/^2, 7:^ - e|- 

L 2m ) 

We set 

{Aer :=9t-0c{0)-mQ,{0). 
Using the triangle inequality, we get the following results: 

||r - 9{t) - {AeriNlll^ = \\a{t) - (Z\0)°°lAr||,oo 

< Wait) - a,{t)\ + \a,{t) - (z40)°°Iw||,oo 

< 2D{a{t)) + \\a,{t) - (Z\0)°°ljv||£-, 



which implies 

This completes the proof. 



lim \\e{t) - {9' + {Ae)°°lN)\\i^ =0. 



3.2 Numerical simulations 

In order to confirm our analytic proof and its validity, we numerically test various situations, which 
is delineated as follows. We employ the 4th order Runge-Kutta method to numerically integrate the 
modified Kuramoto model, Eq. where 



the time step At = 10"^ and N = 100. 
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Fig. 1 (Color online) Semi-logarithmic plots of the relative error, against time, t, which show 

the exponential decay for various conditions: A'^ = 100 and 100 realizations from {Oi(0)} C [— yfj, jq], and 
{cji(O)} C [— TT, tt]. We test two different distributions of natural frequencies, (a) for the uniform distribution. 
Framework A (left) with (m=0.1, K=2.0) and Framework B (right) with (m=0.2, K=2.0) and (b) for 
Cauchy (Lorentzian) distrilsution, Framework A (left) with (m=0.004, 7^=40.0) and Framework B (right) 
with (m=0.01, _ft:=40.0). 



Table 2 The values in the center of mass frame for selected three initial conditions among 100 realizations 
for Figs. [2 and El 



Sample No. 


Oc(0) 




1 


0.020411279668 


-0.056797709228 


2 


0.015377250725 


0.283280670045 


3 


0.022909592396 


-0.004959481400 



As the natural frequency distribution, we consider 

gc{Q) = : Cauchy distribution, gu{^) ~ — '■ uniform distribution. 

One Half of natural frequencies are randomly chosen from the given analytic distribution form in 
the positive part J7 > 0, and the other half is set to be negative of chosen values, so that the mean 
of f2i set to be 0. Initial configurations of {{0^(0)}, {^^(0) = ^,(0)}} are uniformly chosen from the 

intervals: 9i € ~ fs' 15 ^^'^ ^ ["'"'' time t = 0, respectively. 

The following procedure is how we obtain numerical results in this paper: Prepare several sets 
of natural frequencies and initial conditions, {{f2i},{9i{0)},{uji{0)}}. Among them, pick up a set, 
||l2iockcd|^ {6)f )(0)}, {ujf\0)}} with which oscillators converge to a phase-locked state. In this pa- 
per, we test various initial setups, including the simplest ones: the case for initially static, a;i(0) = 
for any oscillators, and the case for Kuramoto-type dynamic, uji{0) = fit + J^j^i sin(0j — 9i). 
We find that for random initial configurations, there are no significant difference between numerical 
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(a) 1 





Re Z; 



(b) 1 





Re Z; 



Fig. 2 (Color online) Three different phase-locked states (o, x, +) on complex plane, which coincide 
with one another by the rotation with the relative phase-shift for A*' — 100, {di{0)} C [— jf;, yj], and 
{ciJi(O)} C [— 7r,7r]. Two different distributions of natural frequencies are tested with the same setups of 
Fig. [Tl summarized in Table [1] The dashed line on the complex plane is guided for eyes as drawing an unit 
circle whose center is located at the origin point. 



results in the steady state. As a result, we can detect phase-locked states as checking the condition 
\9i{t + dt) — 9i{t)\ < £ at every time step, where we set e = 10""'^^. If all phases satisfy the given 
condition, we assume that the system reaches the neighborhood of phase-locked states. 

In Fig. [1] we first check out the 4th statement in Remark 2.2 with Framework A and Framework 

B as stated in Sec. 3, where without loss of generality, we take 6{Q) G ~ fgi xfj a-nd a;(0) G [— tt, tt] 

and measure the relative error from one phase-locked state to one another as follows: 

£{t) := £{t) - f,(0), £{t) 0{t) - 0{t), f,(0) (^^(0) - ^c(O) + J7i(cj,(0) - ci,(0)))lw, 

where X and correspond to the perturbed value from one selected setting and the value in 
the "center of mass" frame, respectively. The values used in our calculation are shown in Table [2] 
As expected, the relative error ||£(t)|| decays exponentially fast to zero, which is the numerical 
verification of Theorem 12. II as well. 

Once {fli] = {fi^"^^^'^} is found, we then fix it as the natural frequency distribution and 
change the initial conditions of {{61,(0)}, {w,(0)}}. With the set {{^2]°'='^°'^}, {6lf ^(0)}, {wf ^(0)}}, 
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Fig. 3 (Color online) For phase-locked states, the analytic estimate (Theorem 1 2. l|l of the relative phase-shift 
is confirmed in Framework A and Framework B with two distributions of natural frequencies: the uniform 
distribution and the Cauchy distribution, respectively. We plot all the cases of Figs. [1] and [21 summarized 
in Table [1] Our numerical simulations and analytical estimates are in perfect agreement with the guide line 
of y = X. 



||l2iockcd| |g,(2)^Q)| |^(2)^Q^||^ and so on. The gathered 100 phase-locked states of final phases, 

{0[^^^}, {ol'^^^},..., {O^^^^^'^}, are plotted as Zi = ^^^1+^^ on the complex plane as shown in Fig. [5] 
where only three cases are taken. Figure [5] shows that different initial conditions reach indeed the 
same phase-locked state with the relative phase-shift, {A6)°° , between one and another, which are 
tested in Framework A and Framework B for two different distributions, respectively. 

The relative phase-shift is the amount of a rotated angle, which is needed to collapse numerical 
data of several phase-locked states with different initial conditions but the same natural frequency 
distribution. Figure |3] shows that all 100 sets of final phases coincide with one another. The sum 
of errors, 

E.!Ii - S^^^l is less than IQ-i", where i is an oscillator index and j, k are sample 
indices. In Fig. [31 numerically estimated values of are compared with analytically estimated 

values (shift prediction) using initial conditions, based on Theorem 12. 1[ which implies that the 
analytic estimate is correct. 

Finally, in Figs. 31 [SI and [H we numerically show that arbitrary initial conditions (neither 
Framework A nor Framework B) also converge to an unique and stable phase- locked state as long 
as keeping the same natural frequency distribution. This implies that the restrictions of Framework 
A and Framework B are not necessary for the system to converge to a phase-locked state for a given 
natural frequency distribution. All the conditions, parameters, values used in our simulation tests 
are shown in Table [T] and Table [2l 



12 




Fig. 4 (Color online) Semi-logarithmic plots of the relative error, against time, t, which show 

the exponential decay for various conditions; A'^ = 100 and 100 realizations from {Oi(0)} C [—jq, jq], and 
{uji{0)} C [— TT, tt]. We test two different distributions of natural frequencies, (a) for the uniform distribution, 
(m=1.0, K=S.O) (left) and (m=0.2, K^2.0) (right), and (b) for the Cauchy distribution, (m=0.5, 7^=40.0) 
(left) and (m=0.1, K=20.0) (right). Note that these parameter setups satisfy neither Framework A nor 
Framework B. 



4 Conclusion 

In summary, we have presented a simple proof for the nonlinear stability of some classes of phase- 
locked states in terms of the modified Kuramoto model with finite inertia in the context of £°°-noTm. 
Phase-locked states of the modified Kuramoto model with f2c — correspond to the equilibrium 
solutions of the system 

K ^ 

1=1 

The solvability of the above nonlinear equation is not clear at all. In the earlier work j20] . 
the existence of phase-locked states are established via a time- asymptotic approach, i.e., instead 
of solving the above nonlinear system (j4.1[) directly, the time-dependent system (|l.ip was solved 
from some admissible initial configurations. Therefore, in the time-asymptotic limit, the desired 
nonlinear system (|4.ip can be solved. This time-asymptotic approach can reveal the fine structure of 
the phase- locked states (see [21] for the original Kuramoto model). Moreover, two frameworks were 
presented for the well-selected parameters and initial configurations to guarantee the validity of this 
time-asymptotic approach. 

In this paper, we have adopted the same framework ideas employed in the aforementioned work, 
where the existence of phase-locked states were only investigated, and showed that the phase-locked 
states whose existence is guaranteed by [20] are orbitally ^°°-stable in the sense that its small 
perturbation leads to the relative phase-shift from the original phase-locked state. Based on this 
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Fig. 5 (Color online) Three different phase-locked states (o, x, +) on complex plane, which coincide 
with one another by the rotation with the relative phase-shift for A*' — 100, {di{0)} C [— jf;, yj], and 
{ciJi(O)} C [— 7r,7r]. Two different distributions of natural frequencies are tested with the same setups of 
Fig. |4l summarized in Table [1] The dashed line on the complex plane is guided for eyes as drawing an unit 
circle whose center is located at the origin point. 



fact, we claim that the phase-locked state has in some sense a robust structure. This implies why we 
can often observe a typical phase-locked states in biological systems. As can be seen in numerical 
simulation tests (see Fig. [2]), Framework A and Framework B in Sec.[2]do not generate all the possible 
phase-locked states and our stability theory do not cover the phase-locked states whose phase- 
diameter is larger than ^. At present, we cannot provide a complete classification for the stability 
of phase-locked states to the system (|l.ip . such as the clear answers of the following questions: For 
a given m,,K and {^2,}, how many phase-locked states exist? Are all phase-locked states orbitally 
stable? Of course, it is impossible to find unstable phase-locked states in our numerical investigation 
because intrinsic numerical errors make unstable phase-locked states invisible. Thus, we leave these 
intriguing issues for future works to be investigated |22| as well as the nature of phase transitions in 
the modified Kuramoto model 

Finally, based on our extensive numerical simulation tests, we also note that the restriction of 
the phase-diameter in phase-locked states is no longer necessary to support the orbital stability, 
D{9^) < ^, which is related to the framework setup. In other words, we numerically observe that 
the phase- locked state with D{6'^) > ^ satisfies its uniqueness and orbital stability as well. 
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Fig. 6 (Color online) For phase-locked states, the analytic estimate (Theorem 1 2. l|l of the relative phase-shift 
is confirmed in Framework A and Framework B with two distributions of natural frequencies: the uniform 
distribution and the Cauchy distribution, respectively. We plot all the cases of Figs. [4] and [5l summarized 
in Table [1] Our numerical simulations and analytical estimates are in perfect agreement with the guide line 
of y = X. 
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